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APPLICATION OF RICHARDSON* S EXTRAPOLATION TO NUMERICAL 


EVALUATION OF SONIC-BOOM INTEGRALS 

By William B. Igoe 
Langley Research Center 


SUMMARY 


An application of Richardson T s extrapolation to the numerical evaluation of 
the sonic-boom integrals occurring in the theory of Whitham (Communications on 
Pure and Applied Mathematics, August 1952) has been considered. Equations have 
been developed for evaluating the sonic-boom integrals, by using second deriva- 
tives which are obtained from the central difference formulas of finite differ- 
ence techniques. It has been shown by error analysis that these equations have 
an order of error proportional to the square of the solution point interval 
spacing for the far-field sonic-boom integral, and to the three-halves power of 
the interval spacing for the near-field sonic-boom integral. Extrapolation of 
the numerical results to zero interval spacing on the basis of the predicted 
order of error is demonstrated, and some numerical examples confirm the essen- 
tial validity of this procedure. Some of the numerical results, however, show 
that the extrapolation technique must be used with caution. 


INTRODUCTION 


The purpose of this report is to present an application of Richardson *s 
extrapolation as a method of improving the accuracy of a numerical evaluation of 
the sonic-boom integrals. These integrals are encountered in using the theory 
of reference 1 for the calculation of the sonic-boom overpressure fields due to 
the flow pattern of a supersonic aircraft or projectile. In general, numerical 
methods must be used in the calculation since the integrals are a function of 
the effective area distribution which is usually known only numerically. The 
lift as well as the volume contribution to the overpressure field of an aircraft 
configuration may be included in the sonic-boom theory through the use of an 
equivalent body concept. Supersonic-area- rule cuts may be used to determine the 
contribution of the aircraft volume to the effective area distribution of the 
equivalent body; thus, the effect of Mach number on these contributions is 
introduced. These extensions of the basic theory of reference 1, as applied to 
the calculation of sonic-boom overpressure fields, have been summarized in 
reference 2. 

The procedure referred to as Richardson's extrapolation was first presented 
in reference and an expanded version was later presented in reference 4. As 
stated in reference 4, the philosophy of this process is to reverse the order of 



passing to the limit in the solution of a problem. In an analytical solution 
involving an application of the methods of the calculus, the passage to the 
limit is taken at the earliest possible stage in the formulation of the problem. 
In a numerical solution, however, this procedure is unworkable, and the solution 
must therefore be obtained in finite differences. In applying Richardson T s 
extrapolation, the approach to the limit is deferred until after the problem has 
been solved for a finite number of discrete increments. If the variation of the 
solution error with increment size is known, the results may then be extrapo- 
lated to what they would be for vanishingly small increments. In references 3 
and 4, the theory of the deferred approach to the limit is applied only to prob- 
lems which are solved by central difference techniques, resulting in what is 

called "h^ extrapolation. " However, the philosophy itself is not necessarily 
restricted to problems of this class. It is the general theory of the deferred 
approach to the limit which is applied in this report, and the particular form 
of the extrapolation used is derived from an analysis of the discretization 
error. 

A number of numerical methods of evaluating the sonic-boom integrals are 
available. (See, for example, refs. 2 and 5 to 7- ) The method of reference 2 
involves the use of the second derivatives of the effective area distribution of 
an equivalent body in a straightforward application of the theory of reference 1 
for both the near- and far- field solutions. Reference 5 presents a method of 
evaluating the near- field overpressure integral involving only first derivatives. 
The far-field overpressure may be readily evaluated once the near- field solution 
is obtained. References 6 and 7 present methods of evaluating the far-field 
overpressure integrals where again only first derivatives are required. The 
numerical procedure for evaluation of the sonic-boom integrals which is utilized 
in this report is essentially the same as that of reference 2, the main differ- 
ence being in the method by which the derivatives of the equivalent body effec- 
tive area distribution are obtained. For this step, the standard central dif- 
ference formulas are used. This method was chosen because it appeared most 
readily amenable to an error analysis and to the subsequent application of 
Richardson’s extrapolation. 

In the following sections of this report, after a description of the symbol 
notation which is used, the near- and far-field overpressure expressions are 
introduced from reference 1 to show the significance of the sonic-boom inte- 
grals. Next the three-point formulas of the standard central difference equa- 
tions are used to obtain the required second derivatives of the equivalent body 
effective area distribution. The second derivatives are then integrated twice 
to obtain a check of the original area distribution. Following this the second 
derivatives are used to obtain expressions for the sonic-boom integrals in a 
form which is suitable for an error analysis. The presentation of the foregoing 
material is essentially a preface to the analysis of the discretization error of 
the given sonic-boom integral solution and the subsequent application of 
Richardson’s extrapolation. The error is obtained as a function of the solution 
interval spacing and only the lowest order terms are retained. Finally some 
numerical examples are presented to show a possible application of the extrapo- 
lation technique. 


2 



SYMBOLS 


All length dimensions are for a body of unit length. 

A(x) effective cross-sectional area of equivalent body which combines 

effects of volume and lift, nondimensionalized with respect to 
square of equivalent body length 

aj coefficient in jth order term of polynomial used to represent A(x) 


p ( x) = i r di 

2n ^0 Vx - | 


h 


width of interval in 


x between regularly spaced body stations 


I(x) = 



d£ 


i integer designating ith derivative 

J integer representing maximum order of polynomial 

j integer representing the order of polynomial terms 

K reflection factor (equals unity if pressure wave does not contact 

reflective surface) 

k arbitrary exponent 

M Mach number 

m,n integers designating body stations starting from zero at nose 

N integer representing maximum station number at base or rear extremity 

of equivalent body 

p reference pressure for uniform atmosphere 

Q quantity representing either F or I 

R^(x),R2(x) radii, nondimensionalized with respect to equivalent body length 

r perpendicular distance from flight path of equivalent body to field 

point at which the sonic-boom overpressure is evaluated, nondimen- 
sionalized with respect to equivalent body length 
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x,| longitudinal distance from nose along equivalent body axis, positive 

downstream, nondimensionalized with respect to equivalent body 
length 

x Q value of x for which integral l(x) is maximum 

0 = )ImP - 1 

y ratio of specific heats, 1.4 for air 

incremental pressure (sonic-boom overpressure) due to flow field of 
equivalent body at supersonic speeds 

& = x - Xn 


Subscripts: 

check first or second integral of numerical approximation of A"(x) 

error difference between exact value and numerical approximation of a 

function 

exact value of function evaluated analytically 

extrap value of function obtained by Richardson's extrapolation 

far field asymptotic value of £p 

m,n,N function evaluated at station m,n,N 

max maximum positive value of a function 

Primes, Roman numerals, and (i) are used as superscripts to indicate deriv- 
atives with respect to x. 


SONIC-BOOM OVERPRESSURE EQUATIONS 


At a sufficiently large distance from the body the incremental pressure at 
a specific point in the flow field of a projectile at supersonic speeds is 
given, to a first approximation, by the theory of reference 1 as 


Ap _ K/M 2 F(x) 
P \|20r 
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where 


F(x) 


1_ 

2 « 



d£ 


The configuration of the shock waves in the supersonic flow field and the magni- 
tude of the pressure increments across the shock waves are determined from area 
balancing requirements which must he applied to the F(x) function. A complete 
description of the influence of the F(x) function with regard to the strength 
and location of the shock waves is given in reference 1. A review and applica- 
tion of this theory is given in reference 8, and a numerical method of accom- 
plishing the area balancing technique as applied to the F(x) function to 
obtain pressure signatures is presented in reference 9 * Under far-field condi- 
tions, reference 1 shows that the pressure increment due to the front shock 
asymptotically approaches the value 



far field 


y(2p) 1 AK[l(x 0 )] 1 / 2 

\jy + lr^A 


( 2 ) 


where l( x 0 ) 


is the maximum positive value of the integral 


I(x) = 



d£ 


The functions F(x) and l(x) represent the sonic-boom integrals which are to 
be evaluated. It can be seen from equations (l) and (2) that when the values of 
the functions F(x) and l(x) are known, the incremental pressures in the flow 
field may then be determined, both under asymptotic and nonasymptotic condi- 
tions. As mentioned previously, reference 2 has presented a numerical method 
(involving second derivatives) for evaluating the functions F(x) and l(x). 

An alternative method for evaluating these functions is developed in a form con- 
sidered suitable for the subsequent analysis of the discretization error. 


NUMERICAL SECOND DERIVATIVE OF AREA DISTRIBUTION 


The second derivatives of the effective area distribution can be approxi- 
mated numerically by using the three-point formula of the standard central dif- 
ference equations 


< “ ft | n ■ ' “» + ^i) 


( 3 ) 
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Ill 


which is equivalent to taking the second derivative of a parabolic arc which is 
passed through three adjacent points. This expression is used to represent the 
value of the second derivative over an interval h in width centered about the 
nth station or in the range of x 

( n ‘ |)h < x < (n + |)h 

At the nose, the second derivative is obtained by 

^“^2 ( oSxS !) W 

Axr 0 h d ' ' 

Equation (4) is strictly applicable only to configurations with pointed noses. 
However, this is not necessarily a restrictive condition to the numerical method 
developed herein since the theory of reference 1 to which it is applied is also 
similarly restricted. At the base of the body, the second derivative in the 
interval at the base is assumed equal to the second derivative in the interval 
immediately preceding it 

%“ a K- 1 ((l.O-!)<*<l.o) (5) 

where the normalized configuration length is unity. This assumption at the base 
is equivalent to applying the forward difference formula of finite difference 
techniques at the base instead of the central difference formula which is 
applied elsewhere on the configuration. A check of the A"(x) approximation is 
obtained by integrating twice to obtain A(x) check . The first integral for the 

check function is 

A ’<*W- / 0 Xa " <5) ai (6) 

which may be integrated with the approximations of equations (3), ( 4 ) , and (5) 
to yield 

A ’Wcheck” V (° SxS t) (7) 

The derivation of equation (8) is shown in detail since the same method is used 
in the summations for A(x) c - [ieck , F(x), and l(x). This method yields expres- 
sions in a form which are well suited to the analysis for the discretization 
error. Equation (6) may be rewritten 
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nh/2 n( n + 2) h px 

A ’( X Wck - / A "^ ^ + LJ, 1N A "^> d * + J, lN A " (l) d& 

u n=l ( n_ 2/ ~ ^ 


( n .i)h 


The second derivative approximations of equations ( 3 ), (4), and ( 5 ) are used to 
obtain 


n=i ( n -2) h (“‘sr 


h/2 


m- 


* rK)\„ 


An d| 


The intervals over which the integrals are taken are changed as follows 


A’(x) 


check 


r *s«- r 1N ^ 

j 0 Jh/2 Si \ (”-|) h 

- f as «\ ♦ r < 

H> I V§> 


dt 


d| 


Integration and rearrangement of terms yields the following form: 


m 


A '^ X ) check *** A 0 X + ^ (Ai Ai-l) x “ ( n 2 ) h ] 

n=l 


((" - |) h S * S (” + |) h ) 

(8) 


The second integral 


A( *> check * f 0 A ’ (l > 


check 


d| 


(9) 


yields 


A(x) 


check 


\ A^x 2 (° = X = |) (10) 
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and 



Values at station or half-station points are obtained by letting x = mh or 
x = ^m + i^h, respectively, in equation (ll). The foregoing check is performed 

to determine how closely the area distribution for which the numerical solution 
is obtained matches the area distribution for which the solution is desired. 


NUMERICAL INTEGRATION FOR F(x) AND l(x) 


The integral for the function F(x) for area distributions possessing con- 
tinuous first derivatives is 

F(x) = i T X (12) 

2lt J 0 \jx - l 

Now using the same stepwise procedure in equation (12) as was used for the area 
check function 

(o<xj|) (13) 


and 
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and at a half-station point, let x = ^m + ijh 


F -i « 
m+= 

2 


>fh 


m 


A o\ m + i + X W ■ A “- 1 ) 'J m ■ n + 1 


n=l 


(16) 


The second integration for the function l(x) is defined as follows 


I(x) = r F(i) 
'J n 


d| 


(17) 


The integration for l(x) is accomplished with the same stepwise procedure as 
before 


i < x > ~ ^ }/z 


(o< xi|) ( 18 ) 


and 


m 


iVl5/2\ 


I(x) ~ h r ° x5/2 + X K - ^-i) [f - ( n - |) h ] 

1 n=l 


(( m - |) h ^ x < (m + |)h) (19) 


Equation (19) for l(x) may be evaluated at a station point by substituting 
x = mh 




h 3/2 


3 Tt 


and at a half- station point by substituting x 


A>3/ 2 - £ K - - n * f) 5/2 

n=l 

. (- * §)h 


( 20 ) 


2 h 


3/2 


m+i 3 * 

2 


A^(m + ^ (An - A^_x)(m - n + l )^ 2 

n=l 


( 21 ) 
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ERROR ANALYSIS 


Among the important sources of error in the numerical solution is the dis- 
cretization error due to the finite-size increments that are used to obtain 
solutions. This error occurs because the second derivative of the effective 
area distribution is approximated numerically and is assumed constant over a 
discrete station interval. At any given station point, the second derivative 

error is of order h 2 but in the ±h/2 interval about the station point, for 
which the second derivative is considered constant, the order of error degrades 

from order h 2 of the central difference formulas to order h of the forward 
difference formulas. The object of the error analysis is to find the lowest 
order of h on which the discretization error of the F and I functions 
depends. All higher orders of h are considered negligible in Richardson's 
extrapolation as h is made vanishingly small. 


The order of error of the second derivative may be shown as follows. The 
station areas A n _]_, A n , and A n+ ^_ are expanded in a Taylor's series centered 

at the arbitrary point x which lies within the interval x n -^<x^x n +^. 

Letting S = x - Xq, there is obtained 


A n _ 1 = A(x) - (h + 5)A'(x) + (h + 5) 2 £1 M _ ( h + 6 )3 

Of ^ ! 


+ 


(h + S) 4 



(h + 5) 5 



+ . . • 


An = A(x) - 5A ’ (x) + 5 2 


5 3 Al»(x) + A ^x) _ s5 AV(x)_ 

31 in 5 : 


and 


A n+1 = A ( x ) + ( h - 5)A'(x) + (h - 5) 2 - ^ + (h - 5)^ A '" | x ^ 

21 31 


+ (h - 8)^ AlV ^ x ^ + (h - S) 5 - V ( x ). + 

4: 5: 


The above expressions are combined in the central difference formula (eq. (3)) 
to obtain 
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a 2 a 

Ax 2 


n ^l 21 L 


h + 5) 2 - 2S 2 + (h 


+ (h - 8)5] + Ai yfe) . |(h + 5)^ 

i 


which may be simplified to 

A 2 A 
Ax 2 


n 


h 2 + 6(x 

= A"(x) - (x - x n ^A ,n (x) + 


8) 2 ] + A^El [_ (h 

28 1 *- + (h - 8)^] + . 


+ S) 5 + 25 5 

. . o(h5)J 



( 22 ) 


In equation (22), the term A"(x) is the true second derivative and the rest 
of the expression is the error. Equation (22) may be rewritten as 


A 2 A 

Ax 2 


n 


A"(x) + A"(x) 


error 


where 


A" ( x ) error = -( x - x n )A'" (x) + + ^ A IV (x) + . . . o(h5) (23) 

The discretization error in the F and I functions may then he written as 

- ox A" (E) 

terror = £ / ~ ^ 

2rt J 0 \Jx - | 

and 


I(x) 


error 



(25) 


In order to proceed with the discretization error analysis it is necessary 
to make some assumption regarding the shape or form of the A(x) curve. A 
polynomial representation of A(x) is used to illustrate the order of error 
of the numerical solution. Use of a ^th degree polynomial retains error terms 

of order h^ in equation (22) and eliminates terms of order h^ and higher. 
This assumption is considered adequate for the purposes of the error analysis. 
The general polynomial can he written 
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«■> • E =£ 

3=o 


and the ith derivative is 


A^(x) = £ 


J- i 


a.tX') - -*- 

(d - i) '• 


For J = k, the polynomial expression becomes 


A(x) - + —■ + a^x + aQ 

T • 3 * 2 m 


( 26 ) 


and the third and fourth derivatives are 


A'" (x) = a^x + a^ 


(27) 


and 


A IV (x) - 


( 28 ) 


Equations (27 ) , ( 28 ), and (23) may he substituted into equation (24) to obtain 

F(x) ^ a*x^/^ + — h^xr*-/^ - — a lx5/2 

x error jr ) 3 3 12 15 ^ 


I [ a 5 h + ( n " |) a 4 h2 ] [ x - ( n " i) 11 ] 172 ) 

n=l J 


^m - |)h < x < (m + ( 29 ) 


Equation (29) may be evaluated at a station point by letting x = mh 
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F_ « - (a,]* 1 / 2 

m, error ^ Y*5 


m 


-!" 3/2 + I (--“*5 


) 1/2 


11=1 


♦ s^/2 


A ^5/2 + i m l/2 + £ (n - i)(» - n ♦ |) 1/2 

n=l 


( 30 ) 


and at a half- station point "by letting x = + i^h 


N) 


error 




+ a^h 5 / 2 


m 


- i (- + i) 3/2 * I (m - n + 1,1/2 

n=l 

Jl ( m + if /2 + i ( m + i ) l/2 

15 \ 2/ 12 \ 2/ 


X (» - §)<- - - - i> i/2 


n=l 


(3D 


To evaluate the error expression for the I function, equation (29) may be 
used in equation (25) to obtain 


I(x) 


error 


' > ^ { I x5/2 + h !_ [ x ■ ( n - iitf 1 ' } 


m 


+ a., \ x 

4 1 35 


i\J5/2l 


7/2 ♦ fi ^ ** l (■-§)[*- M>] 

n=l 


(( m - |) h = X = ( m + J) 11 ) ^ 
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When I(x) is evaluated at a station point x = mh, equation (32) becomes 

error 


error ~ £ ^ S a 3 h ' 


, 1/2 


- | 1,5/2 + I “ + if ' 2 

n=l 


+ a 4 h^ 


. 1 m 7/2 + _L m 3/2 + £ ( n . i)( m . n + i^3/2 

n=l 


(33) 


and at a half- station point, where x = ^m + equation (32) becomes 


('■*!). 


i- h 2 (a,!, 1 / 2 

3a ) 3 


error 


m 


- | (m + + Yj ( m “ n + 1 ) 


3/2 


n=l 


+ a^/ 2 


35 


(-•if *s<-*r 


^ (n - |)(m - n + l) 5 / 2 
n=l ' 


(3^) 


If approximations for the sums of the series are introduced in equa- 
tions (30), (31) j (33 ) > and (3*0 then, for m large, on an order of magnitude 
basis, the F and I error expressions are reduced to the following 
expressions: 


* 


m, error 


in it ( 

2k * V 


a 3 + 


a 4 x ) 


(35) 


( 





3 h3/2 

2b a 



+ a 4 x ) 


T « 

^m, error 


xV 2 h 2 

2kn 




(36) 

(37) 
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( a 3 “ 2a^x) 


(58) 



xV g h g 

2kn 


Equations (35) and (37) show that the F and I functions, as numerically 
evaluated at a station point using equations (15) and (20) , have discretization 

errors essentially of order h^/^ and i? , respectively. By reducing h, the 
station interval width (and increasing N, the number of station points in a 
solution), the discretization error can be reduced but only at the expense of 
increased computing time and accumulated error. Another method of reducing the 
discretization error is to use Richardson's extrapolation to extrapolate the 
computed values to vanishingly small values of h if the order of error of the 
solution is known. 

To use Richardson's extrapolation, the F or the I functions are eval- 
uated at common station points for two different values of h (h-^ and hg>) • 

Then extrapolated values may be computed by 


p = Fm,! 13 ^ 2 ~ F m,2 h l^ 2 

m, extrap h^^ 2 ^ 3/2 


t = A 1 ^ 2 - T m ,2^1 2 

m, extrap ^2 _ h] 2 


(59) 




Use of this extrapolation should give an improved approximation to the correct 
value provided that round-off and other accumulated errors are small and that 
the station intervals are sufficiently small so as to minimize the effect of the 
higher orders of h which have been neglected in the derivation of the error 
terms. The half-station error results (eqs. (36) and (58)) are of little inter- 
est with regard to extrapolation since it is usually difficult to obtain solu- 
tions at very many matching stations. It may be observed here that the 

A(x) c h eC ]c equation (eq. (ll)) also has a discretization error of order h?. 


STEPS IN NUMERICAL SOLUTION 


A possible procedure for obtaining numerical solutions for the sonic-boom 
equations may be as follows: 

(l) Determine the effective equivalent body area distribution including 
both the volume and lift effects following the methods of reference 2. 
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(2) Obtain values at the required station points selecting an interval 

spacing which will give the smallest desired value of h. Select multiples of 
the smallest interval spacing which will give common station points for larger 
values of h (at least two spacings h^ and h2 where h2 = 2h-^, for 

instance). 

(5) Compute numerical approximations of the second derivative A"(x) by- 
using equations (3), (4), and (5). 

(4) Check the resulting equi valent area distribution A(x) chec j c Uy using 
equations (10) and (ll). 

(5) Compute the I’m and I m values at all required station points for all 
desired station interval spacings by using equations (15) and (20). 

(6) At common station points for the various station interval spacings 
selected, apply the extrapolation formulas of equations (39) and (4o) to the 

I’m and Im values, respectively, or plot the results against h 3/2 

and h^, 

respectively, and extrapolate graphically to h = 0. In the event that computed 
results are obtained which do not follow the order of error shown here, deter- 
mine the order of error h k by graphical means, if possible, and complete the 
extrapolation graphically or use the formula 


1, extrap 


- Sn,2 h l k 

hg k - hx k 


(41) 


where Q is a quantity representing either the F or the I values. 

(7) Use methods such as those presented in references 8 or 9 to compute the 
desired pressure fields. 

(8) If the maximum far-field overpressure is desired from equation (2), use 
graphical or analytical interpolation to determine l( x ) max * One simple analyt- 
ical interpolation scheme is to fit a parabola through the three adjacent larg- 
est positive computed values of 1^. If these three consecutive values are des- 
ignated Ii, I2 , and Ij, then l(x) Tnfly is given approximately by 


l(x) « I-. 
max 1 


1 (-? 1 ! + jgg - Ilf 

8 (I X - 212 + I 5 ) 


(42) 


In general, better results are obtained from equation (42) if it is applied to 
extrapolated results at the necessarily wider spacing rather than to the origi- 
nal computed results at the closest spacing. The computed values of l( x ) max 
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cannot be extrapolated very well since they will not necessarily be at common 
station points or x values. 


NUMERICAL EXAMPLES 


Parabolic Bodies 

Two numerical examples for a parabolic body are shown. The first is for a 
simple smooth body with a parabolic distribution of radii 

RiCx) = 0. 7725 |l - 4(x - 0.5) 2 ] (43) 

The second is for a body similar to the first one hut with a hump equivalent to 
adding the amount of cross-sectional area contained in a body with parabolic 
distribution of radii 


^(x) = 0.252 [l - 100 (x - 0.4) 2 ] (0.3 < x < O.5) (44) 

Only normal or cross-sectional areas (corresponding to M = 1.0 cuts) were used 
for the equivalent body area distribution in the numerical examples. The 
resulting normal area distributions are for the smooth parabolic body 


and for the parabolic body 

AgU) 

Ag(x) 

-^(x) 


'x 

OJH 

P5 

II 

'x 

< 


with bump 


= rtR^(x) 

“X 

(o < x < O.j) 

= jiJr|(x) + R|(x)J 

(0.3 < x < 0.5) 

'x 

f 

II 

(0.5 < X < 1.0) 


(45) 


(46) 


The enclosed volume of the parabolic body with bump is approximately 2 percent 
greater than the volume of the smooth parabolic body. These bodies were chosen 
only for illustrative purposes and do not necessarily resemble any realistic 
configuration. It should be observed that the smooth parabolic body area dis- 
tribution function (eq.. (45)) is a 4th degree polynomial in x and therefore 
corresponds to what was assumed in the error analysis. The results of the 
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error analysis (eqs. (35) to ( 38 )) are therefore most nearly applicable to the 
solution for this body. The solution for the parabolic body with bump however 
will contain higher order terms. For the smooth parabolic body, the constants 
in equation ( 2 6 ) are: 


Sq = 0.0 

a-]_ = 0.0 

&2 = 19. 0962 jr 
a 3 = -114.5772* 
a^ = 229-1544* 


The A check , F m , and 3^ values, computed with equations (ll), (15) .> 

and (20), are compared with the exact values in figure 1 for both bodies for an 
interval spacing of h = 0.1 and h = 0.05. It may be seen from figure 1 that 
increasing the number of station points from 10 to 20 improved the solution 
accuracy. A further illustration of the improvement in solution accuracy with 
number of station points (N = l/h) is shown in figure 2. Here, F m and 3^ 

are shown in an extrapolation against h 5/2 and h?, respectively, at a body 

station point of x = 0.4 in figure 2(a) and at x = 0.8 in figure 2(b) for 

both bodies. All but the F function for the smooth parabolic body at x = 0.4 

show an orderly progression. In this one plot, the F scale is greatly 
expanded and even the N = 10 (h = 0.1) solution is quite accurate (error less 

than 1 percent). Here it is possible that errors of order other than h^/^ or 
sources of error other than the discretization error are influencing the behav- 
ior of the solutions. 


Logarithmic plots of the absolute value of the error fractions 


F, 


m, error 


exact 


and 


m, error 


obtained by subtracting the exact values from the approximate 


■*-exact 

confuted values, are presented in figure 3 as a function of h for selected 
body station points. These data show, by the nearly constant slope of most of 
the lines, that the order of error is approximately as predicted by equa- 
tions ( 35 ) and ( 37 ) • A few notable exceptions are for x = 0.3 and x = 0.5 
for the parabolic body with bump, and for x = 1.0 for both bodies. The first 
two stations correspond to the points where the bump starts and ends, and 
although the first derivative of the area distribution is continuous here, the 
second derivative is not. The same is true at the base (x = 1.0) of both 
bodies. At such points, apparently, the order of error is not predictable. 
However, graphically at the bump discontinuities, the order of error seems to 
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approach h ' for the F values and hr' for the I values. Error frac- 
tions for I I are presented in figure 4 for data computed by using the half- 
m+i 

station formula of equation (21) and verify the h 2 order of error in the half- 
station error expression of equation (38). 

Computed values of the I function at body stations of x = O.32, 
x = 0.3 4, and x = O.36 for h intervals of 0.01 and 0.02 are shown in the 
following table for the smooth parabolic body: 


h 

Ija (smooth parabolic body) for - 


x = 0.32 

x = 0.54 

x = O.36 

0.02 

0.85643415 

O.86271865 

0.86095690 

.01 

.85781819 

.86407731 

.86228505 


These values were used in the extrapolation formula of equation ( 40 ) and the 
extrapolated results are also shown as follows: 


x = 0.32 

0 . 8582795 ^ 


The value I max = 0.86483733 at x Q = 0.3455236 was obtained by interpolation 
with equation ( 42 ) and the values of I e xtrap This compares with an 

exact value of ^max, exact = 0.86483645 at x = 0.3^5^915 with an error in 
Imax one part in a million. 


I , for - 

extrap 

x = 0.34 

x = 0. 36 

0.86453020 

0.86272777 


Practical Airplane Shape 

Figure 5 shows a comparison of the F m and I m functions computed by the 

central differences method used in this report (circular symbols) and by the 
method of reference 2 (square symbols). The area distribution chosen was one 
for a practical airplane shape similar to the one presented in figure 2 of ref- 
erence 2 . The comparison shows that the two methods, both of which involve 
integrals containing the second derivative explicitly, predict substantially 
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the same results for "both the F and I functions. In the comparison of 
Acheck* a ^- so in figure the square symbols correspond to the input values 

since a feature of the method of reference 2 is that, at the station points, 
the input area is recovered exactly when the second derivative is integrated 
twice for the area. However, departures from the desired area distribution 
would occur in between the station points because of the oscillatory nature of 
the second derivatives in that method. 

The convergence of the solutions with decreasing interval spacing is shown 
in figure 6 for body stations x = 0.4 and x = 0.8. Extrapolation of the 

I function against h^ appears quite regular but extrapolation of the 

F function against apparently indicates that convergence is not obtained 

until the interval spacing is quite small (h = 0.01, h = 0.02). Convergence 
would not be expected until the interval spacing is sufficiently small so as to 
define the waviness or bumpiness of the area distribution correctly. Probably 
a minimum of four stations within a given wave or bump would be required to give 
adequate definition and only solutions with this interval spacing and closer 
would be expected to show convergence in an extrapolation. It is therefore 
apparent that the extrapolation technique must be used with caution, and cannot 
be applied indiscriminately or without some prior knowledge of the nature of 
convergence of a solution with decreasing interval spacing. 


CONCLUDING REMARKS 


An application of Richardson’s extrapolation to the numerical evaluation 
of the sonic-boom integrals occurring in the theory of Whitham (Communications 
on Pure and Applied Mathematics, August 1952) has been considered. Equations 
have been developed for evaluating the sonic-boom integrals, using second deriv- 
atives which are obtained from the central difference formulas of finite differ- 
ence techniques. It has been shown by error analysis that these equations have 
an order of error proportional to the square of the solution point interval 
spacing for the far-field sonic-boom integral, and to the three-halves power of 
the interval spacing for the near-field sonic-boom integral. Extrapolation of 
the numerical results to zero interval spacing on the basis of the predicted 
order of error is demonstrated, and some numerical examples confirm the essen- 
tial validity of the procedure. Some of the numerical results however show 
that the extrapolation technique must be used with caution. 


Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va. , September 6, 1966, 
126-16-03-04-25. 
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Figure 5.- Comparison of numerical values of A chec ( < , F m , and l m for a practical airplane shape computed by the central difference method 

and by the method of reference 2 for h = 0.01. 
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